Contrasting resistance and resilience to light variation of the coupled oxic and anoxic components of an experimental microbial ecosystem
Marcel Suleiman
7/27/2021
Introduction
This markdown contains all analyses done for the manuscript “Contrasting resistance and resilience of the coupled oxic and anoxic components of an experimental microbial community”. The first part deals with the oxygen data, the second part deals with the microbial 16S rRNA amplicons (PacBio full length sequencing).
Short summary of the project:
-8 micro-ecosystems, incubated at 24 °C for 8 days under light:dark cycles of 16:8h -afterwards, 4 micro-ecosystems (stressor treatment) were incubated in darkness for 7 days -afterwards, the stressor-treated columns were incubated again as the controls -Oxygen was measured at the top and bottom of the liquid part continuously every 5 min -samples were taken at day 8 (prior-stressor treatment), day 15 ( stressor-treatment), day 19 (short time recovery) and day 35 (long term recovery) at the height of the top and bottom oxygen sensor, respectively, for 16S rRNA full length amplicon sequencing
Aim of the study: Analysis of resistance and resilience of the oxygen concentration and the oxic/anoxic microbial community
Oxygen
data handling
oxygen_data <- read_excel("oxygen_data")Fix the date and times
oxygen_data <- oxygen_data %>%
mutate(Date_time = mdy_hms(paste0(Date, Time))) %>%
select(-Date, -Time)Calculate hourly average
hourly_value <- oxygen_data %>%
mutate(Date = date(Date_time),
Hour = hour(Date_time)) %>%
group_by(Sensor_Name, Date, Hour) %>%
dplyr::summarise(mean_Value = mean(Value)) %>%
mutate(Date_time = ymd_h(paste0(Date, Hour)))## `summarise()` has grouped output by 'Sensor_Name', 'Date'. You can override using the `.groups` argument.## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.
## Warning: All formats failed to parse. No formats found.hourly_value <- hourly_value %>% mutate(Day = as.numeric(Date - ymd("2020-12-02")))Calculate the DAILY mean, minimum, and maximum of the hourly means
daily_values <- hourly_value %>%
ungroup() %>%
group_by(Date, Sensor_Name) %>%
dplyr::summarise(mean_O2 = mean(mean_Value),
min_O2 = min(mean_Value),
max_O2 = max(mean_Value),
amplitude_O2 = max_O2 - min_O2) %>%
mutate(Date_time_midpoint = ymd_hms(paste0(Date, "12:00:00"))) %>%
pivot_longer(names_to = "Variable", values_to = "Oxygen", cols = 3:6)## `summarise()` has grouped output by 'Date'. You can override using the `.groups` argument.put in missing values for the days where there is no data
missing_data <- crossing(Sensor_Name = unique(daily_values$Sensor_Name),
Variable = unique(daily_values$Variable),
Date = dmy(paste(c(24:27), "12-2020")),
Date_time_midpoint = ymd_hms(paste0(Date, " 12:00:00")))%>%
mutate(Oxygen = NA)
daily_values <- rbind(daily_values, missing_data) %>%
mutate(Sensor_Name=tolower(Sensor_Name)) %>%
separate(Sensor_Name, into= c("Treatment", "Sensorposition","Column"),
remove=FALSE)
daily_values <- daily_values %>% mutate(Day = as.numeric(Date - ymd("2020-12-02")),
Variable2 = factor(Variable, levels = c("mean_O2", "max_O2", "min_O2", "amplitude_O2")))plot the data
hourly values
plot_oxygen_a <- hourly_value %>%
mutate(Day2 = Day*24 + Hour) %>%
ggplot(aes(x = Day2, y = mean_Value, col = Sensor_Name)) +
geom_line() +
geom_rect(xmin = 195, xmax = 360,ymin = 0, ymax = 100, col = "lightblue", fill= "lightblue") +
geom_line() +
geom_rect(xmin = 510, xmax = 640,ymin = 0, ymax = 100, col = "white", fill= "white") +
scale_color_manual(values =c("Control_top_1" = "black", "Control_top_2" = "black", "Control_top_3" = "black", "Control_top_4" = "black","Disturbance_top_1" ="red","Disturbance_top_2" ="red","Disturbance_top_3" ="red","Disturbance_top_4" ="red","Control_bottom_1" = "grey", "Control_bottom_2" = "grey", "Control_bottom_3" = "grey", "Control_bottom_4" = "grey","Disturbance_bottom_1" ="grey","disturbance_bottom_2" ="grey","disturbance_bottom_3" ="grey","Disturbance_bottom_4" ="grey")) +
theme_bw() + theme(panel.grid = element_blank()) + scale_x_continuous(breaks = seq(0, 817, by=48), labels=seq(0,35, by=2))+
theme(legend.position = "none") +
ylab("Hourly mean oxygen value [%]") +
xlab("Day") +
geom_point(data=data.frame(Day2=c(190, 360, 456, 840), mean_Value=rep(-1, 4)), colour="black", size=3, shape=4)
plot(plot_oxygen_a)
just top sensor with all 4 variables and daily mean
plot_oxygen_b<- daily_values %>%
filter(Sensor_Name %in% c("control_top_1","control_top_2","control_top_3","control_top_4","disturbance_top_1","disturbance_top_2","disturbance_top_3","disturbance_top_4")) %>%
ggplot(aes(x = Day, y = Oxygen, col = Sensor_Name)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_line() +
facet_wrap( ~ Variable2, ncol = 4) +
scale_color_manual(values =c("control_top_1" = "black", "control_top_2" = "black", "control_top_3" = "black", "control_top_4" = "black","disturbance_top_1" ="red","disturbance_top_2" ="red","disturbance_top_3" ="red","disturbance_top_4" ="red","control_bottom_1" = "black", "control_bottom_2" = "black", "control_bottom_3" = "black", "control_bottom_4" = "black","disturbance_bottom_1" ="red","disturbance_bottom_2" ="red","disturbance_bottom_3" ="red","disturbance_bottom_4" ="red")) +
theme_bw() + theme(panel.grid = element_blank())+
theme(legend.position = "none") +
geom_point(data=data.frame(Day=c(7.5, 15, 19, 35), Oxygen=rep(-7, 4)), shape=4 ,colour="black", size=2) +
ylab("Daily mean oxygen value [%]") +
xlab("Day")
plot(plot_oxygen_b)
Applying t-test for every day to get p.value for top sensors
analysis <- daily_values %>%
filter(Sensor_Name != "control_top_1") %>%
na.omit() %>%
group_by(Date, Variable, Date_time_midpoint, Sensorposition) %>%
do(m1 = tidy(lm(Oxygen ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest() %>%
filter(term=="Treatmentdisturbance", Sensorposition == "top")## Warning: `cols` is now required when using unnest().
## Please use `cols = c(m1)`analysis <- analysis %>% mutate(Day = as.numeric(Date - ymd("2020-12-02")),
Variable2 = factor(Variable, levels = c("mean_O2", "max_O2", "min_O2", "amplitude_O2")))Plot the estimates
#final
plot_estimate_oxygen<- analysis %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_ribbon(aes(ymin=conf.low,ymax=conf.high, x=Day), alpha=.4)+
geom_rect(xmin = 22, xmax = 25,
ymin = -Inf, ymax = 100, col = "white", fill="white", alpha=1) +
facet_wrap( ~ Variable2, ncol = 4) +
scale_color_manual(values =c("control_top_1" = "black", "control_top_2" = "black", "control_top_3" = "black", "control_top_4" = "black","disturbance_top_1" ="red","disturbance_top_2" ="red","disturbance_top_3" ="red","disturbance_top_4" ="red","control_bottom_1" = "black", "control_bottom_2" = "black", "control_bottom_3" = "black", "control_bottom_4" = "black","disturbance_bottom_1" ="red","disturbance_bottom_2" ="red","disturbance_bottom_3" ="red","disturbance_bottom_4" ="red")) +
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
geom_point(data=data.frame(Day=c(7.5, 15, 19, 35), estimate=rep(-25, 4)), colour="black", size=2, shape=4) +
ylab("Estimate") +
xlab("Day")
print(plot_estimate_oxygen)
Microbial community
Load phyloseq data, remove chloroplast sequences, just include samples from Experiment part III
ps <- readRDS("ps_oxic_anoxic.rds")
ps_new <- subset_taxa(ps, Order != "Chloroplast")
taxa_names(ps_new) <- paste0("Seq", seq(ntaxa(ps)))Calculation of alpha diversity index
diversity_test<- estimate_richness(ps_new, split=TRUE, measures=NULL)## Warning in estimate_richness(ps_new, split = TRUE, measures = NULL): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.diversity_test$ID <- rownames(diversity_test)
names <- read_excel("barcodes_FKIII.xlsx")## New names:
## * `` -> ...1names1 <- as.data.frame(names)
names1$ID <- gsub("--", "..", names1$...1)
diversity_analysis <- merge(diversity_test, names1, by=c("ID"))
stderr <- function(x, na.rm=FALSE) {
if (na.rm) x <- na.omit(x)
sqrt(var(x)/length(x))
}
div_means <- diversity_analysis %>%
dplyr::filter(Position %in% c("bottom", "top")) %>%
dplyr::group_by(Position, Day, Treatment) %>%
dplyr::summarize(Shannon_mean = mean(Shannon),
Shannon_SE = stderr(Shannon))## `summarise()` has grouped output by 'Position', 'Day'. You can override using the `.groups` argument.div_raw <- diversity_analysis
plot_a_diversity_top <- ggplot() +
geom_rect(data=div_raw, aes(xmin = 8, xmax = 15,
ymin = -Inf, ymax = Inf), col = "lightblue", fill="lightblue", alpha=0.1) +
geom_point(data = subset(div_raw, Position == "top"), aes(x=Day, y=Shannon, col=Treatment,shape=Position), size=3, alpha=.7)+
scale_shape_manual(values=c(17,16))+
geom_point(data=subset(div_means, Position == "top"), aes(x=Day, y=Shannon_mean, col=Treatment), shape="+", size = 10, position = "dodge", width = 0.25) +
geom_line(data=subset(div_means, Position =="top"), aes(x=Day, y=Shannon_mean, col=Treatment), linetype="dashed") +
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
theme_bw() + theme(panel.grid = element_blank())+
theme(legend.position = "none") +
scale_x_continuous(limits=c(0,36),breaks = c(0,8,15,19,35)) +
ggtitle("top layer community") +
theme(plot.title = element_text(size = 11))## Warning: Ignoring unknown parameters: widthprint(plot_a_diversity_top)## Warning: Width not defined. Set with `position_dodge(width = ?)`
plot_a_diversity_bottom <- ggplot() +
geom_rect(data=div_raw, aes(xmin = 8, xmax = 15,
ymin = -Inf, ymax = Inf), col = "lightblue", fill="lightblue", alpha=0.1) +
geom_point(data = subset(div_raw, Position == "bottom"), aes(x=Day, y=Shannon, col=Treatment,shape=Position), size=3, alpha=.7)+
scale_shape_manual(values=c(17,16))+
geom_point(data=subset(div_means, Position == "bottom"), aes(x=Day, y=Shannon_mean, col=Treatment), shape="+", size = 10, position = "dodge", width = 0.25) +
geom_line(data=subset(div_means, Position =="bottom"), aes(x=Day, y=Shannon_mean, col=Treatment), linetype="dashed") +
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
theme_bw() + theme(panel.grid = element_blank())+
theme(legend.position = "none") +
scale_x_continuous(limits=c(0,36),breaks = c(0,8,15,19,35)) +
ggtitle("bottom layer community")+
theme(plot.title = element_text(size = 11))## Warning: Ignoring unknown parameters: widthprint(plot_a_diversity_bottom)## Warning: Width not defined. Set with `position_dodge(width = ?)`
estimates of the alpha diversity analysis
analysis_shannon <- div_raw %>%
group_by(Day, Position) %>%
do(m1 = tidy(lm(Shannon ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
# plot the estimates
plot_estimates_shannon_top <- analysis_shannon%>%
filter(Position=="top") %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(ymin=conf.low,ymax=conf.high, x=Day),width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
ylab("Estimate") +
xlab("Day") +
scale_x_continuous(limits=c(0,37),breaks = c(0, 8,15,19,35))
print(plot_estimates_shannon_top)
plot_estimates_shannon_bottom <- analysis_shannon%>%
filter(Position=="bottom") %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(ymin=conf.low,ymax=conf.high, x=Day),width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
ylab("Estimate") +
xlab("Day") +
scale_x_continuous(limits=c(0,37),breaks = c(0, 8,15,19,35))
print(plot_estimates_shannon_bottom)
#Anova
## top community
diversity_analysis_top<- diversity_analysis %>%
filter(Position=="top")
lm_treatment <- lm(Shannon ~ Treatment*Day, data = diversity_analysis_top)
autoplot(lm_treatment)## Warning: `arrange_()` was deprecated in dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
anova(lm_treatment)## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 0.0138 0.0138 0.0423 0.8385271
## Day 1 6.5483 6.5483 20.0874 0.0001225 ***
## Treatment:Day 1 1.6860 1.6860 5.1719 0.0311192 *
## Residuals 27 8.8018 0.3260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(lm_treatment)##
## Call:
## lm(formula = Shannon ~ Treatment * Day, data = diversity_analysis_top)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5282 -0.3754 0.1258 0.4857 0.6453
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.51565 0.31188 17.685 2.25e-16 ***
## TreatmentDisturbance -0.84603 0.44627 -1.896 0.0687 .
## Day -0.06889 0.01441 -4.782 5.47e-05 ***
## TreatmentDisturbance:Day 0.04647 0.02044 2.274 0.0311 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.571 on 27 degrees of freedom
## Multiple R-squared: 0.4838, Adjusted R-squared: 0.4264
## F-statistic: 8.434 on 3 and 27 DF, p-value: 0.0004077##top sensor shannon index is affected "postponed", and it does not show resilience
Anova_treatment <- lm(Shannon ~ Treatment*as.factor(Day), data = diversity_analysis_top)
autoplot(Anova_treatment)
anova(Anova_treatment)## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 0.0138 0.0138 0.0988 0.756165
## as.factor(Day) 3 11.7919 3.9306 28.1259 7.133e-08 ***
## Treatment:as.factor(Day) 3 2.0299 0.6766 4.8418 0.009353 **
## Residuals 23 3.2143 0.1398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(Anova_treatment)##
## Call:
## lm(formula = Shannon ~ Treatment * as.factor(Day), data = diversity_analysis_top)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.68415 -0.19311 0.01744 0.30636 0.44183
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3096 0.1869 28.406 < 2e-16 ***
## TreatmentDisturbance -0.4739 0.2643 -1.793 0.08616 .
## as.factor(Day)15 -0.6177 0.2643 -2.337 0.02852 *
## as.factor(Day)19 -1.9468 0.2643 -7.365 1.72e-07 ***
## as.factor(Day)35 -1.9154 0.2643 -7.246 2.25e-07 ***
## TreatmentDisturbance:as.factor(Day)15 0.0423 0.3891 0.109 0.91436
## TreatmentDisturbance:as.factor(Day)19 0.8408 0.3738 2.249 0.03438 *
## TreatmentDisturbance:as.factor(Day)35 1.1886 0.3738 3.179 0.00418 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3738 on 23 degrees of freedom
## Multiple R-squared: 0.8115, Adjusted R-squared: 0.7541
## F-statistic: 14.14 on 7 and 23 DF, p-value: 5.529e-07AIC(lm_treatment, Anova_treatment)## df AIC
## lm_treatment 5 58.94428
## Anova_treatment 9 35.71636#bottom community
diversity_analysis_bottom<- diversity_analysis %>%
filter(Position=="bottom")
##bottom sensor shannon index is affected directly, and it does show resilience
lm_treatment <- lm(Shannon ~ Treatment*Day, data = diversity_analysis_bottom)
autoplot(lm_treatment)
anova(lm_treatment)## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 6.0279 6.0279 7.9775 0.008631 **
## Day 1 22.0637 22.0637 29.1999 9.209e-06 ***
## Treatment:Day 1 0.2105 0.2105 0.2786 0.601757
## Residuals 28 21.1570 0.7556
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(lm_treatment)##
## Call:
## lm(formula = Shannon ~ Treatment * Day, data = diversity_analysis_bottom)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7771 -0.6599 0.0257 0.7631 1.4147
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.85624 0.47482 8.121 7.67e-09 ***
## TreatmentDisturbance 1.18319 0.67150 1.762 0.08898 .
## Day -0.07561 0.02193 -3.448 0.00181 **
## TreatmentDisturbance:Day -0.01637 0.03102 -0.528 0.60176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8693 on 28 degrees of freedom
## Multiple R-squared: 0.5722, Adjusted R-squared: 0.5264
## F-statistic: 12.49 on 3 and 28 DF, p-value: 2.307e-05#bottom sensor shannon index is affected directly, and it does show resilience
Anova_treatment <- lm(Shannon ~ Treatment*as.factor(Day), data = diversity_analysis_bottom)
autoplot(Anova_treatment)
anova(Anova_treatment)## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 6.0279 6.0279 12.050 0.001978 **
## as.factor(Day) 3 25.2738 8.4246 16.841 4.204e-06 ***
## Treatment:as.factor(Day) 3 6.1515 2.0505 4.099 0.017537 *
## Residuals 24 12.0059 0.5002
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(Anova_treatment)##
## Call:
## lm(formula = Shannon ~ Treatment * as.factor(Day), data = diversity_analysis_bottom)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.22400 -0.31305 0.00521 0.39579 1.33089
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.39590 0.35364 9.603 1.08e-09 ***
## TreatmentDisturbance -0.06541 0.50012 -0.131 0.897026
## as.factor(Day)15 -0.73628 0.50012 -1.472 0.153962
## as.factor(Day)19 -1.14224 0.50012 -2.284 0.031520 *
## as.factor(Day)35 -2.10233 0.50012 -4.204 0.000315 ***
## TreatmentDisturbance:as.factor(Day)15 1.83580 0.70728 2.596 0.015862 *
## TreatmentDisturbance:as.factor(Day)19 1.78258 0.70728 2.520 0.018783 *
## TreatmentDisturbance:as.factor(Day)35 0.11542 0.70728 0.163 0.871739
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7073 on 24 degrees of freedom
## Multiple R-squared: 0.7573, Adjusted R-squared: 0.6865
## F-statistic: 10.7 on 7 and 24 DF, p-value: 4.674e-06AIC(lm_treatment, Anova_treatment)## df AIC
## lm_treatment 5 87.57164
## Anova_treatment 9 77.44136Calculation rel abundance, just sequences that appear > 0.1 %, making dataframes for each column position
ps_rel <- transform_sample_counts(ps_new, function(x) x / sum(x) )
relab_threshold <- 0.001
ps_relab <- filter_taxa(ps_rel, function(x) !(sum(x < relab_threshold) == length(x)), TRUE)
ntaxa(ps_new)## [1] 9028ntaxa(ps_relab)## [1] 1856ps_relative <- transform_sample_counts(ps_relab, function(x) x / sum(x))
df_rel <- psmelt(ps_relative)
df_rel_top <- df_rel %>%
filter(Position == "top")
df_rel_bottom <- df_rel %>%
filter(Position == "bottom")NMDS
NMDS Plot
mds_whole <- ps_relative@otu_table %>%
as.data.frame() %>%
metaMDS(.,
distance = "bray",
k = 2, ## number of dimensions to reduce to
try = 200, ## number of random starts to try
autotransform = FALSE ## best not to use
)## Run 0 stress 0.1530756
## Run 1 stress 0.1530827
## ... Procrustes: rmse 0.00148772 max resid 0.0103001
## Run 2 stress 0.1530831
## ... Procrustes: rmse 0.001531195 max resid 0.01032903
## Run 3 stress 0.1571977
## Run 4 stress 0.236323
## Run 5 stress 0.1571977
## Run 6 stress 0.1571977
## Run 7 stress 0.1530827
## ... Procrustes: rmse 0.001485717 max resid 0.01028305
## Run 8 stress 0.1530829
## ... Procrustes: rmse 0.001513705 max resid 0.01031943
## Run 9 stress 0.1530827
## ... Procrustes: rmse 0.001483084 max resid 0.01026715
## Run 10 stress 0.1530827
## ... Procrustes: rmse 0.00148431 max resid 0.01027275
## Run 11 stress 0.2078723
## Run 12 stress 0.1530827
## ... Procrustes: rmse 0.001483782 max resid 0.01027015
## Run 13 stress 0.153083
## ... Procrustes: rmse 0.001527468 max resid 0.01032705
## Run 14 stress 0.153083
## ... Procrustes: rmse 0.001522597 max resid 0.0103238
## Run 15 stress 0.1530827
## ... Procrustes: rmse 0.001484558 max resid 0.01027117
## Run 16 stress 0.1530827
## ... Procrustes: rmse 0.001486136 max resid 0.01027742
## Run 17 stress 0.1571977
## Run 18 stress 0.1530827
## ... Procrustes: rmse 0.001486199 max resid 0.01028388
## Run 19 stress 0.1571977
## Run 20 stress 0.1571814
## Run 21 stress 0.1571977
## Run 22 stress 0.2275183
## Run 23 stress 0.1530827
## ... Procrustes: rmse 0.00148717 max resid 0.01029209
## Run 24 stress 0.1530827
## ... Procrustes: rmse 0.00148708 max resid 0.0102902
## Run 25 stress 0.1835705
## Run 26 stress 0.1530829
## ... Procrustes: rmse 0.001522374 max resid 0.01032379
## Run 27 stress 0.1530756
## ... Procrustes: rmse 2.585429e-05 max resid 0.0001694419
## ... Similar to previous best
## Run 28 stress 0.1530827
## ... Procrustes: rmse 0.001485381 max resid 0.01027804
## Run 29 stress 0.1863566
## Run 30 stress 0.1571814
## Run 31 stress 0.1530756
## ... New best solution
## ... Procrustes: rmse 1.911943e-05 max resid 0.0001243164
## ... Similar to previous best
## Run 32 stress 0.1530756
## ... Procrustes: rmse 4.391557e-05 max resid 0.0002895892
## ... Similar to previous best
## Run 33 stress 0.2462748
## Run 34 stress 0.2452603
## Run 35 stress 0.1530756
## ... Procrustes: rmse 2.010439e-05 max resid 0.0001319844
## ... Similar to previous best
## Run 36 stress 0.1530827
## ... Procrustes: rmse 0.00148587 max resid 0.01028955
## Run 37 stress 0.1530756
## ... Procrustes: rmse 1.716225e-05 max resid 0.0001121669
## ... Similar to previous best
## Run 38 stress 0.1571977
## Run 39 stress 0.227854
## Run 40 stress 0.1571977
## Run 41 stress 0.1530827
## ... Procrustes: rmse 0.001486417 max resid 0.01029439
## Run 42 stress 0.1530756
## ... Procrustes: rmse 1.800349e-05 max resid 0.0001184586
## ... Similar to previous best
## Run 43 stress 0.2437621
## Run 44 stress 0.153083
## ... Procrustes: rmse 0.001520742 max resid 0.01031471
## Run 45 stress 0.1530756
## ... Procrustes: rmse 2.018693e-05 max resid 0.0001329001
## ... Similar to previous best
## Run 46 stress 0.1571977
## Run 47 stress 0.1530756
## ... Procrustes: rmse 2.90842e-06 max resid 1.408669e-05
## ... Similar to previous best
## Run 48 stress 0.1530827
## ... Procrustes: rmse 0.001482603 max resid 0.01026852
## Run 49 stress 0.1530827
## ... Procrustes: rmse 0.001485844 max resid 0.01028957
## Run 50 stress 0.1530827
## ... Procrustes: rmse 0.00148322 max resid 0.01027252
## Run 51 stress 0.1530827
## ... Procrustes: rmse 0.001485107 max resid 0.01028072
## Run 52 stress 0.1530829
## ... Procrustes: rmse 0.001515176 max resid 0.01032225
## Run 53 stress 0.1530831
## ... Procrustes: rmse 0.001532307 max resid 0.0103358
## Run 54 stress 0.153083
## ... Procrustes: rmse 0.001529385 max resid 0.01033284
## Run 55 stress 0.1530827
## ... Procrustes: rmse 0.001482715 max resid 0.01027004
## Run 56 stress 0.183571
## Run 57 stress 0.1571977
## Run 58 stress 0.153083
## ... Procrustes: rmse 0.001501185 max resid 0.01030609
## Run 59 stress 0.1530756
## ... Procrustes: rmse 2.806648e-06 max resid 1.687471e-05
## ... Similar to previous best
## Run 60 stress 0.1530827
## ... Procrustes: rmse 0.001486085 max resid 0.01028057
## Run 61 stress 0.1530831
## ... Procrustes: rmse 0.00154512 max resid 0.01034472
## Run 62 stress 0.1571977
## Run 63 stress 0.1530756
## ... New best solution
## ... Procrustes: rmse 9.745351e-06 max resid 6.388499e-05
## ... Similar to previous best
## Run 64 stress 0.1530827
## ... Procrustes: rmse 0.001483923 max resid 0.01027472
## Run 65 stress 0.2383096
## Run 66 stress 0.1530757
## ... Procrustes: rmse 0.0001074537 max resid 0.0007097788
## ... Similar to previous best
## Run 67 stress 0.2078723
## Run 68 stress 0.153083
## ... Procrustes: rmse 0.001533951 max resid 0.01033842
## Run 69 stress 0.1530828
## ... Procrustes: rmse 0.001502232 max resid 0.01030918
## Run 70 stress 0.1530827
## ... Procrustes: rmse 0.00148633 max resid 0.01028869
## Run 71 stress 0.1530827
## ... Procrustes: rmse 0.001487739 max resid 0.01028381
## Run 72 stress 0.1571977
## Run 73 stress 0.1530756
## ... Procrustes: rmse 1.115265e-05 max resid 7.306532e-05
## ... Similar to previous best
## Run 74 stress 0.1530832
## ... Procrustes: rmse 0.001553023 max resid 0.01035205
## Run 75 stress 0.1530827
## ... Procrustes: rmse 0.001481748 max resid 0.0102606
## Run 76 stress 0.1530828
## ... Procrustes: rmse 0.001498745 max resid 0.01030342
## Run 77 stress 0.1571977
## Run 78 stress 0.1530827
## ... Procrustes: rmse 0.001488154 max resid 0.01029276
## Run 79 stress 0.1530827
## ... Procrustes: rmse 0.001492443 max resid 0.01029632
## Run 80 stress 0.1530827
## ... Procrustes: rmse 0.001488339 max resid 0.01030626
## Run 81 stress 0.1530756
## ... Procrustes: rmse 1.113055e-05 max resid 7.311672e-05
## ... Similar to previous best
## Run 82 stress 0.1530829
## ... Procrustes: rmse 0.001509633 max resid 0.0103157
## Run 83 stress 0.1530756
## ... Procrustes: rmse 7.169728e-06 max resid 4.657897e-05
## ... Similar to previous best
## Run 84 stress 0.1530827
## ... Procrustes: rmse 0.00148912 max resid 0.01028799
## Run 85 stress 0.226658
## Run 86 stress 0.1530828
## ... Procrustes: rmse 0.001508742 max resid 0.01031494
## Run 87 stress 0.1530827
## ... Procrustes: rmse 0.00148184 max resid 0.01026222
## Run 88 stress 0.2288104
## Run 89 stress 0.1530831
## ... Procrustes: rmse 0.001544824 max resid 0.01035264
## Run 90 stress 0.1530831
## ... Procrustes: rmse 0.001545692 max resid 0.01035093
## Run 91 stress 0.1530831
## ... Procrustes: rmse 0.001537873 max resid 0.0103442
## Run 92 stress 0.1530756
## ... Procrustes: rmse 2.447473e-06 max resid 1.156627e-05
## ... Similar to previous best
## Run 93 stress 0.1530827
## ... Procrustes: rmse 0.001483834 max resid 0.01027246
## Run 94 stress 0.1530828
## ... Procrustes: rmse 0.00150334 max resid 0.01031119
## Run 95 stress 0.2383482
## Run 96 stress 0.1530829
## ... Procrustes: rmse 0.00151533 max resid 0.01032149
## Run 97 stress 0.153083
## ... Procrustes: rmse 0.00153258 max resid 0.01033407
## Run 98 stress 0.1530827
## ... Procrustes: rmse 0.001490857 max resid 0.01029209
## Run 99 stress 0.1530829
## ... Procrustes: rmse 0.001515011 max resid 0.01032188
## Run 100 stress 0.2044065
## Run 101 stress 0.2284918
## Run 102 stress 0.1530827
## ... Procrustes: rmse 0.00148422 max resid 0.01027434
## Run 103 stress 0.1835701
## Run 104 stress 0.1530827
## ... Procrustes: rmse 0.001487628 max resid 0.01029498
## Run 105 stress 0.1530829
## ... Procrustes: rmse 0.001519059 max resid 0.01032514
## Run 106 stress 0.1530827
## ... Procrustes: rmse 0.001482414 max resid 0.010266
## Run 107 stress 0.1571977
## Run 108 stress 0.2233197
## Run 109 stress 0.153083
## ... Procrustes: rmse 0.00153533 max resid 0.01034335
## Run 110 stress 0.1530828
## ... Procrustes: rmse 0.001504003 max resid 0.01031102
## Run 111 stress 0.1571977
## Run 112 stress 0.2198539
## Run 113 stress 0.1530827
## ... Procrustes: rmse 0.001487103 max resid 0.01029527
## Run 114 stress 0.1530829
## ... Procrustes: rmse 0.001510092 max resid 0.01031824
## Run 115 stress 0.1530756
## ... Procrustes: rmse 5.588646e-06 max resid 3.533489e-05
## ... Similar to previous best
## Run 116 stress 0.1530757
## ... Procrustes: rmse 8.080733e-05 max resid 0.0005236894
## ... Similar to previous best
## Run 117 stress 0.1530756
## ... Procrustes: rmse 1.849313e-06 max resid 1.105e-05
## ... Similar to previous best
## Run 118 stress 0.1530756
## ... Procrustes: rmse 2.52858e-06 max resid 1.035051e-05
## ... Similar to previous best
## Run 119 stress 0.1571977
## Run 120 stress 0.1530828
## ... Procrustes: rmse 0.001505058 max resid 0.01031287
## Run 121 stress 0.1571977
## Run 122 stress 0.236946
## Run 123 stress 0.2201603
## Run 124 stress 0.1571814
## Run 125 stress 0.1530831
## ... Procrustes: rmse 0.001543113 max resid 0.0103587
## Run 126 stress 0.1530756
## ... Procrustes: rmse 8.588398e-06 max resid 5.180093e-05
## ... Similar to previous best
## Run 127 stress 0.1835683
## Run 128 stress 0.1571977
## Run 129 stress 0.1530756
## ... Procrustes: rmse 1.594459e-06 max resid 4.629072e-06
## ... Similar to previous best
## Run 130 stress 0.1530831
## ... Procrustes: rmse 0.001539571 max resid 0.01036053
## Run 131 stress 0.1530831
## ... Procrustes: rmse 0.001537212 max resid 0.01033847
## Run 132 stress 0.1530827
## ... Procrustes: rmse 0.001483814 max resid 0.01027562
## Run 133 stress 0.2044065
## Run 134 stress 0.1530756
## ... Procrustes: rmse 1.198952e-05 max resid 7.876109e-05
## ... Similar to previous best
## Run 135 stress 0.2078903
## Run 136 stress 0.1864996
## Run 137 stress 0.1571977
## Run 138 stress 0.1530827
## ... Procrustes: rmse 0.001483878 max resid 0.01027624
## Run 139 stress 0.1530756
## ... Procrustes: rmse 1.102587e-05 max resid 7.150153e-05
## ... Similar to previous best
## Run 140 stress 0.2418881
## Run 141 stress 0.2149961
## Run 142 stress 0.2247954
## Run 143 stress 0.1571814
## Run 144 stress 0.1571814
## Run 145 stress 0.1530827
## ... Procrustes: rmse 0.00148572 max resid 0.01028685
## Run 146 stress 0.1571977
## Run 147 stress 0.1530827
## ... Procrustes: rmse 0.00148448 max resid 0.01027579
## Run 148 stress 0.1530829
## ... Procrustes: rmse 0.001512976 max resid 0.01031603
## Run 149 stress 0.1835696
## Run 150 stress 0.1571977
## Run 151 stress 0.2202608
## Run 152 stress 0.1530756
## ... Procrustes: rmse 8.212996e-06 max resid 5.359464e-05
## ... Similar to previous best
## Run 153 stress 0.1530756
## ... Procrustes: rmse 7.548606e-06 max resid 4.950594e-05
## ... Similar to previous best
## Run 154 stress 0.1530827
## ... Procrustes: rmse 0.001486845 max resid 0.01029449
## Run 155 stress 0.1530827
## ... Procrustes: rmse 0.001480886 max resid 0.01024511
## Run 156 stress 0.1530829
## ... Procrustes: rmse 0.001516893 max resid 0.010325
## Run 157 stress 0.1530827
## ... Procrustes: rmse 0.001479907 max resid 0.01024182
## Run 158 stress 0.2203492
## Run 159 stress 0.1530827
## ... Procrustes: rmse 0.001487245 max resid 0.01029409
## Run 160 stress 0.1530827
## ... Procrustes: rmse 0.001485919 max resid 0.01028685
## Run 161 stress 0.1530757
## ... Procrustes: rmse 0.0001006575 max resid 0.0006602985
## ... Similar to previous best
## Run 162 stress 0.1530827
## ... Procrustes: rmse 0.001487802 max resid 0.01029371
## Run 163 stress 0.1530829
## ... Procrustes: rmse 0.00150922 max resid 0.01030585
## Run 164 stress 0.1530828
## ... Procrustes: rmse 0.001495956 max resid 0.01030119
## Run 165 stress 0.2348738
## Run 166 stress 0.153083
## ... Procrustes: rmse 0.001530306 max resid 0.01033334
## Run 167 stress 0.1571977
## Run 168 stress 0.1571977
## Run 169 stress 0.1530831
## ... Procrustes: rmse 0.001540548 max resid 0.01034203
## Run 170 stress 0.1530827
## ... Procrustes: rmse 0.001483741 max resid 0.01027234
## Run 171 stress 0.1530756
## ... Procrustes: rmse 4.863105e-06 max resid 3.150449e-05
## ... Similar to previous best
## Run 172 stress 0.1571977
## Run 173 stress 0.1530756
## ... Procrustes: rmse 1.542525e-06 max resid 7.538724e-06
## ... Similar to previous best
## Run 174 stress 0.2044065
## Run 175 stress 0.153083
## ... Procrustes: rmse 0.001528137 max resid 0.01033225
## Run 176 stress 0.1530829
## ... Procrustes: rmse 0.001509756 max resid 0.01031674
## Run 177 stress 0.1530831
## ... Procrustes: rmse 0.001545012 max resid 0.01034572
## Run 178 stress 0.2427732
## Run 179 stress 0.1530827
## ... Procrustes: rmse 0.001492616 max resid 0.01029683
## Run 180 stress 0.1530827
## ... Procrustes: rmse 0.001483253 max resid 0.01027038
## Run 181 stress 0.2427561
## Run 182 stress 0.1530831
## ... Procrustes: rmse 0.001547099 max resid 0.01034647
## Run 183 stress 0.2198537
## Run 184 stress 0.1571814
## Run 185 stress 0.1530756
## ... Procrustes: rmse 9.468755e-06 max resid 6.206903e-05
## ... Similar to previous best
## Run 186 stress 0.1571977
## Run 187 stress 0.1530756
## ... New best solution
## ... Procrustes: rmse 1.888246e-06 max resid 1.151079e-05
## ... Similar to previous best
## Run 188 stress 0.1530756
## ... Procrustes: rmse 1.28233e-05 max resid 8.444567e-05
## ... Similar to previous best
## Run 189 stress 0.1530757
## ... Procrustes: rmse 5.772538e-05 max resid 0.0003430395
## ... Similar to previous best
## Run 190 stress 0.2277348
## Run 191 stress 0.2450719
## Run 192 stress 0.2220019
## Run 193 stress 0.1530827
## ... Procrustes: rmse 0.001486224 max resid 0.01028027
## Run 194 stress 0.2430565
## Run 195 stress 0.1530828
## ... Procrustes: rmse 0.001507977 max resid 0.0103144
## Run 196 stress 0.153083
## ... Procrustes: rmse 0.001525161 max resid 0.01033065
## Run 197 stress 0.2044065
## Run 198 stress 0.1530827
## ... Procrustes: rmse 0.001486139 max resid 0.01028887
## Run 199 stress 0.2373475
## Run 200 stress 0.1571977
## *** Solution reached## 0.16
mds_whole_res <- ps_relative @sam_data %>%
as.tibble() %>%
select(Treatment, Column, Day, Position,Day_name) %>%
bind_cols(as.tibble(scores(mds_whole, display = "sites")))## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 objectnmds_plot_supplements <- ggplot(mds_whole_res, aes(x = NMDS1, y = NMDS2)) +
geom_point(aes(shape = Position, color=Treatment),
size = 2) +
scale_shape_manual(values =c(17,16))+
facet_wrap("Day", ncol=2)+
theme_bw() + theme(panel.grid = element_blank())+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))
print(nmds_plot_supplements)
NMDS separated for NMDS1 and NMDS2 axis
NMDS_top <- mds_whole_res %>%
filter(Position=="top")
NMDS_bottom <- mds_whole_res %>%
filter(Position=="bottom")
# Dynamics
plot_NMDS1_top <- NMDS_top %>%
ggplot(aes(x = Day, y = NMDS1, col = Treatment)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point(shape=16, size = 2)+
geom_smooth() +
theme_bw() + theme(panel.grid = element_blank())+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
ggtitle("top layer community") +
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
theme(legend.position = "none")+
theme(plot.title = element_text(size = 11))
print(plot_NMDS1_top)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
plot_NMDS2_top <- NMDS_top %>%
ggplot(aes(x = Day, y = NMDS2, col = Treatment)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point(shape=16, size = 2) +
geom_smooth() +
theme_bw() + theme(panel.grid = element_blank())+
ggtitle("top layer community")+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
theme(legend.position = "none")+
theme(plot.title = element_text(size = 11))
print(plot_NMDS2_top)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
plot_NMDS1_bottom <- NMDS_bottom %>%
ggplot(aes(x = Day, y = NMDS1, col = Treatment)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point(shape=17, size = 2) +
geom_smooth() +
theme_bw() + theme(panel.grid = element_blank())+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
ggtitle("bottom layer community")+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
theme(legend.position = "none")+
theme(plot.title = element_text(size = 11))
print(plot_NMDS1_bottom)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
plot_NMDS2_bottom<- NMDS_bottom %>%
ggplot(aes(x = Day, y = NMDS2, col = Treatment)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point(shape=17, size = 2) +
geom_smooth() +
theme_bw() + theme(panel.grid = element_blank())+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
ggtitle("bottom layer community")+
theme(legend.position = "none")+
theme(plot.title = element_text(size = 11))
print(plot_NMDS2_bottom)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
estimate analysis for NMDS data
analysis_NMDS1_top <- NMDS_top %>%
group_by(Day, Day_name) %>%
do(m1 = tidy(lm(NMDS1 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
plot_estimate_NMDS1_top <- analysis_NMDS1_top %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
ylab("Estimate") +
xlab("Day")
print(plot_estimate_NMDS1_top)
analysis_NMDS2_top <- NMDS_top %>%
group_by(Day, Day_name) %>%
do(m1 = tidy(lm(NMDS2 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
plot_estimate_NMDS2_top <- analysis_NMDS2_top %>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
geom_hline(yintercept=0, linetype="dashed") +
ylab("Estimate") +
xlab("Day")
print(plot_estimate_NMDS2_top)
analysis_NMDS1_bottom <- NMDS_bottom %>%
group_by(Day, Day_name) %>%
do(m1 = tidy(lm(NMDS1 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
plot_estimate_NMDS1_bottom <- analysis_NMDS1_bottom%>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed")+
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
ylab("Estimate") +
xlab("Day")
print(plot_estimate_NMDS1_bottom)
analysis_NMDS2_bottom <- NMDS_bottom %>%
group_by(Day, Day_name) %>%
do(m1 = tidy(lm(NMDS2 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
unnest(cols = c(m1)) %>%
filter(term=="TreatmentDisturbance")
plot_estimate_NMDS2_bottom <- analysis_NMDS2_bottom%>%
ggplot(aes(x = Day, y = estimate)) +
geom_rect(xmin = 8, xmax = 15,
ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
geom_point() +
geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
theme_bw() + theme(panel.grid = element_blank())+
geom_hline(yintercept=0, linetype="dashed") +
scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
ylab("Estimate") +
xlab("Day")
print(plot_estimate_NMDS2_bottom)
final plot of figures combined
p1 <- plot_oxygen_a
p2 <- plot_oxygen_b
p3 <- plot_a_diversity_top
p4<- plot_a_diversity_bottom
p5 <- plot_estimates_shannon_top
p6 <- plot_estimates_shannon_bottom
p7 <- plot_NMDS1_bottom
p8 <- plot_NMDS2_bottom
p9 <- plot_NMDS1_top
p10 <- plot_NMDS2_top
p11 <- plot_estimate_oxygen
p12<- plot_estimates_shannon_top
p13<- plot_estimates_shannon_bottom
p14 <- plot_estimate_NMDS1_bottom
p15 <- plot_estimate_NMDS2_bottom
p16 <- plot_estimate_NMDS1_top
p17 <- plot_estimate_NMDS2_top
#Figure 1
plot_1 <- p1/
p2/
p11
plot_1 +
plot_layout(heights=c(2,1,1)) + plot_annotation(tag_levels = "a")
#Figure 2
plot <- (p4|p13)/
(p3|p12)/
(p7|p14)/
(p8|p15)/
(p9 |p16)/
(p10 |p17)
plot + plot_annotation(tag_levels = "a")## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Width not defined. Set with `position_dodge(width = ?)`## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
#Supplementary Information
Visualization of rel.abundance for supplementary information
df_class<- df_rel %>%
group_by(Class, Treatment,Column,Day, Position) %>%
dplyr::summarize(Abundance=sum(Abundance))## `summarise()` has grouped output by 'Class', 'Treatment', 'Column', 'Day'. You can override using the `.groups` argument.index <- which(df_class$Abundance>=0.05)
class_to_keep <- unique(df_class[index,"Class"])
class_to_keep <- unname(unlist(class_to_keep))
df_class$Class_filter <- ifelse(df_class$Class %in% class_to_keep, df_class$Class,"Zother")
cols <-c("Actinobacteria"="bisque1", "Alphaproteobacteria" = "yellow", "Babeliae"="lightcyan1","Bacilli" = "green4", "Bacteroidia" ="maroon1", "Campylobacteria" ="royalblue1","Zother"="black","Chlorobia" ="thistle1", "Clostridia"= "firebrick1", "Cyanobacteriia" = "magenta4", "Fimbriimonadia" = "grey","Gammaproteobacteria"="red", "Microgenomatia"="grey38","Oligoflexia"="snow1", "Phycisphaerae"="cornflowerblue", "Planctomycetes"="lightskyblue1","Spirochaetia"="lavender", "Vampirivibronia" = "lightsalmon", "Verrucomicrobiae" = "slateblue3")
supplement_plot1 <-df_class%>%
ggplot(aes_string(x = "Column", y = "Abundance", fill="Class_filter" )) +
geom_bar(stat = "identity", position = "stack", col="black") +
facet_wrap(Position~Day, ncol = 4) +
theme(axis.text.x = element_text(angle = -45, hjust = 0, vjust=0.5, size=9))+
ylab("Relative abundance") +
scale_fill_manual("legend", values = cols)
print(supplement_plot1)
Correlation last day, shannon vs.oxygen amplitude
new <- div_raw %>%
filter(Day=="35", Position=="top")
new2 <- daily_values %>%
filter(Day=="35", Variable=="amplitude_O2", Sensorposition=="top") %>% mutate(Column2 = paste0(str_to_title(Treatment),"_", Column))
new3 <- merge(new, new2, by.x=c("Column"), by.y=c("Column2"))
ggplot(new3) + geom_point(aes(x=log(Shannon),y=Oxygen, colour=Treatment.x),size=3, cex.lab=10) +
theme_bw() + theme(panel.grid = element_blank()) +
scale_color_manual(values =c("Control" = "black", "Disturbance" = "red")) +
stat_smooth(aes(x=log(Shannon),y=Oxygen), method="lm") +
ylab("Oxygen (Amplitude)") +
xlab("log(Shannon-index)") +
theme(legend.position = "none") +
theme(axis.text.x = element_text(size=16),axis.text.y =element_text(size=20)) +
theme(axis.title = element_text(size=16))## Warning: Ignoring unknown parameters: cex.lab## `geom_smooth()` using formula 'y ~ x'
cor.test(log(new3$Shannon),new3$Oxygen)##
## Pearson's product-moment correlation
##
## data: log(new3$Shannon) and new3$Oxygen
## t = -5.1134, df = 6, p-value = 0.002192
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9822793 -0.5410140
## sample estimates:
## cor
## -0.9018627#without column 3
new3_no_control_3 <- new3 %>%
filter(Column != "Control_3")
ggplot(new3_no_control_3) + geom_point(aes(x=log(Shannon),y=Oxygen, colour=Treatment.x),size=4, cex.lab=10) +
theme_bw() + theme(panel.grid = element_blank()) +
scale_color_manual(values =c("Control" = "black", "Disturbance" = "red")) +
stat_smooth(aes(x=log(Shannon),y=Oxygen), method="lm") +
ylab("Oxygen (Amplitude)") +
xlab("log(Shannon-index)") +
theme(legend.position = "none") +
theme(axis.text.x = element_text(size=16),axis.text.y =element_text(size=20)) +
theme(axis.title = element_text(size=16))## Warning: Ignoring unknown parameters: cex.lab## `geom_smooth()` using formula 'y ~ x'
cor.test(log(new3_no_control_3$Shannon),new3_no_control_3$Oxygen)##
## Pearson's product-moment correlation
##
## data: log(new3_no_control_3$Shannon) and new3_no_control_3$Oxygen
## t = -2.6205, df = 5, p-value = 0.04707
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.96242875 -0.01788358
## sample estimates:
## cor
## -0.7606971nmds <- mds_whole_res %>%
filter(Day=="35", Position=="top")
nmds2 <- merge(nmds, new2, by.x=c("Column"), by.y=c("Column2"))
ggplot(nmds2, aes(x=NMDS1,y=Oxygen, colour=Treatment.x)) + geom_point()
cor.test(nmds2$NMDS1,nmds2$Oxygen)##
## Pearson's product-moment correlation
##
## data: nmds2$NMDS1 and nmds2$Oxygen
## t = 0.45373, df = 6, p-value = 0.666
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5994779 0.7859367
## sample estimates:
## cor
## 0.1821357ggplot(nmds2, aes(x=NMDS2,y=Oxygen, colour=Treatment.x)) + geom_point()
cor.test(nmds2$NMDS2,nmds2$Oxygen)##
## Pearson's product-moment correlation
##
## data: nmds2$NMDS2 and nmds2$Oxygen
## t = 1.1638, df = 6, p-value = 0.2887
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3949661 0.8705568
## sample estimates:
## cor
## 0.4291481